#!/usr/bin/perl -w
use strict;

# for projecting alignments from protein to R/DNA space
use Bio::Align::Utilities qw(aa_to_dna_aln);
# for input of the sequence data
use Bio::SeqIO;
use Bio::AlignIO;
use Bio::DB::Fasta;

use Getopt::Long;

my $aln_format = 'fasta';

GetOptions(
	   'f|format:s' => \$aln_format,
	   );

# when you run this, do 
# perl protein_to_cds_alignment.pl family.aln cds 

# if you want to use a different alignment format you run it like this
# perl protein_to_cds_alignment.pl -f clustalw family.aln cds
# OR
# perl protein_to_cds_alignment.pl --format clustalw family.aln cds


my ($alignment_file, $cds_database) = @ARGV;

# let's process this
# read in the alignment
my $dbh = Bio::DB::Fasta->new($cds_database);

my $in = Bio::AlignIO->new(-format => $aln_format,
			   -file   => $alignment_file);
if( my $aa_aln = $in->next_aln ) {
    my %seqs;
    for my $s ( $aa_aln->each_seq ) {
	my $id = $s->id;
	if( $id =~ /^tmes/) { # skip tmes
	    $aa_aln->remove_seq($s);
	    next;
	}
	$seqs{$id} = $dbh->get_Seq_by_acc($id);
    }
    my $out = Bio::AlignIO->new(-format => $aln_format,
				-file   => ">$alignment_file.cdsaln");
    my $dna_aln = aa_to_dna_aln($aa_aln, \%seqs);
    
    $out->write_aln($dna_aln);

} else {
    warn("Cannot read an alignment from file '$alignment_file'\n");
}
     
